// GriewankFunction.cpp: implementation of the GriewankFunction class.
//
//////////////////////////////////////////////////////////////////////
#include <iomanip>
#include <fstream>
#include <string>
#include "griewank_function.h"
#include "DataHolder.h"



//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

GriewankFunction::GriewankFunction(int i_NInput)
{
  NInput = i_NInput;
  nCall = 0;
  nMult = 1;

  newA(mCosArray, NInput);          // stores the cosine values
  newA(mSinArray, NInput);          // stores the sine values
  newA(dDerivative, NInput);        // a dummy derivative used in constructing the derivative
  newA(mDerivative,NInput);

  Result.Size(NInput);
}


GriewankFunction::~GriewankFunction()
{
  deleteA(mCosArray);
  deleteA(mSinArray);
  deleteA(dDerivative); 
  deleteA(mDerivative);
}

// takes a point and if it is beyong boundary move the point
// into the boundary box
void GriewankFunction::checkBoundary(vector<double> &Input)
{
  for (unsigned int i = 0; i < Input.size(); i++)
  {
    // reflective boundary
    /*while( fabs(Input[i]) >= 100.0 )
    {
      if ( Input[i] >= 100.0 )
      {
        Input[i] = 100.0 - fabs(Input[i]-100.0);
      }
      else if ( Input[i] <= -100.0 )
      {
        Input[i] = -100.0 + fabs(-100.0-Input[i]);
      }
    }*/

    // periodic boundary
    if ( Input[i] >= 100.0 )
    {
      Input[i] = -100.0 + fmod(Input[i],100.0);
    }
    else if ( Input[i] <= -100.0 )
    {
      Input[i] = 100.0 + fmod(Input[i],100.0);
    }
  }
}

// takes a point and if it is beyong boundary move the point
// into the boundary box
void GriewankFunction::checkBoundary(double *Input)
{
  for (int i = 1; i <= NInput; i++)
  {
    // periodic boundary
    if ( Input[i] >= 100.0 )
    {
      Input[i] = -100.0 + fmod(Input[i],100.0);
    }
    else if ( Input[i] <= -100.0 )
    {
      Input[i] = 100.0 + fmod(Input[i],100.0);
    }
  }
}

//////////////////////////////////////////////////////////////////////
// functions
//////////////////////////////////////////////////////////////////////

double GriewankFunction::Value(double *Input)
{
  // first check that point is not beyong any boundaries
  checkBoundary(Input);

  // Initialise vars answer and danswer (both local)
  double answer = 0.0;
  double danswer = 1.0;

  // Makes answer equal to it self added to
  // the sum of x^2/4000 for each dimension
  for (int m = 1; m <=NInput; m++)
  {
    answer += (Input[m]*Input[m])/4000;
  }

  // Makes danswer equal to itself multiplied
  for(int n=1; n<=NInput; n++)
  {
    danswer *= cos(Input[n]/sqrt(double (n)));
  }

  // sets answer to be equal to itself - danswer + 1 as required
  answer -= danswer;
  answer += 1;

  Result.Input(answer, Input);

  /*if (nCall == (nMult*10))
  {
    cout << "After " << nCall << " function calls, the func value is " << answer << endl;
    nMult++;
  }*/

  nCall++;

  return (answer); 
}

// calculates/builds the function value
//
// f(x) = sum_i=1^n (x_i^2/4000) - prod_i=1^n cos(x_i/sqrt(i)) + 1
//
double GriewankFunction::Value(vector<double> &Input)
{
  // first check that point is not beyong any boundaries
  checkBoundary(Input);

  // Initialise 
  double answer = 0.0;   // to cal sum term
  double danswer = 1.0;  // to cal product term

  // Constructing first term of equation
  for (int m = 0; m < NInput; m++)
  {
    answer += (Input[m]*Input[m])/4000;
  }

  // calculate second product term
  for(int n=0; n < NInput; n++)
  {
    danswer *= cos(Input[n]/sqrt(double (n+1)));
  }

  answer -= danswer; // substract product term
  answer += 1; 

  Result.Input(answer, Input);

  nCall++;

  return (answer); 
}

// calculates/builds the derivative
//
// f(x) = sum_i=1^n (x_i^2/4000) - prod_i=1^n cos(x_i/sqrt(i)) + 1
//
void GriewankFunction::BuildDerivative(double *Input, double *Derivative)
{
  // NInput = dimension of function
  double mCos = 1.0;

  // Constructing equation
  for(int i=1; i<=NInput; i++)
  {
    Derivative[i] = 0.0;
    mCosArray[i] = cos(Input[i]/sqrt(double (i)));
  }

  for(int m=1; m<=NInput; m++)	  
  {
    dDerivative[m] = (Input[m]/2000);  
    mSinArray[m] = (1.0/sqrt(double (m)))*sin(Input[m]/sqrt(double (m)));

    // this is pretty pointless as due to rounding error - mCosArray never reaches zero just becomes
    // very close - however doesn't hurt to be safe.
    if(mCosArray[m] == 0.0)
    {
      cout << "Setting mCosArray[" << m << "] = 1.0";
      mCosArray[m] = 1.0;
    }

    mCos *= mCosArray[m]; 
  }

  for(int n=1; n<=NInput; n++)
  {
    // just reiterating the same point as before however doesn't hurt to be safe.
    if(mCosArray[n] == 0.0)
    {
      cout << "oops - dividing by zero";
    }

    // stores the derivative as an array
    // first element = derivative with respect to the first var
    Derivative[n] = dDerivative[n] + (mSinArray[n]*mCos)/mCosArray[n];
  }

  for(int mw=1; mw<=NInput; mw++)
  {
    mDerivative[mw] = Derivative[mw];
  }
}

// used to test if the function is working properly
void GriewankFunction::test_function()
{
  ofstream out("output/test_of_griewank.txt");
  
  out << "NUMBER_OF_DIMENSIONS: " << NInput << endl 
    << "      Subscript[x,k]      |     Value    |         Derivative[k]" << endl << endl;

  int NRun = 3;
  double *mInput;
  double *mDeriv;
  newA(mInput, NInput);
  newA(mDeriv, NInput);

  for(int run=1; run<=NRun; run++)
  {
    for(int n=1; n<=NInput; n++)
    {
      mInput[n]=(200.0*ran1())-100.0;
    }

    for(int k=1; k<=NInput; k++)
    {
      out << setw(11) << mInput[k];
    }
      // d[1] == 0, doesn't give derivative = 0
      // as d[1] is never 0.
      out << setw(16) << Value(mInput);

      BuildDerivative(mInput, mDeriv);

    for(int l=1; l<=NInput; l++)
    {
      out << setw(16) << mDerivative[l];
    }

    out << "\n\n";
    
  }

  deleteA(mInput);
  deleteA(mDeriv);
  out.close();
}

// used to print the solution (this is the old version, the new version is below)
void GriewankFunction::print_sol(double *solution, double funcVal, char filename[])
{
  ofstream out(filename);	

  out << "Number of dimensions: " << NInput << endl;
  out << "Final Griewank value = " << funcVal << endl;
  out << "Solution coordinates:" << endl;

  for (int i = 1; i <= NInput; i++)
  {
    out << solution[i] << "  ";
  }
 
  out << endl;

  out.close();
}

// used to print the solution to a file
void GriewankFunction::print_sol(char filename[])
{
  ofstream out(filename);	

  // getBestValue() returns the best value stored in the tracker
  out << "Number of dimensions: " << NInput << endl;
  out << "Final Griewank value = " << Result.getBestValue() << endl;
  out << "Solution coordinates:" << endl;

  // prints out the solution coordinates to a file
  Result.printSolution(out);

  out.close();
}


// distance from know global minimum which is located at the origin
double GriewankFunction::distanceFromGlobalMinimum(double *solution)
{
  double toReturn = 0.0;

  for (int i = 1; i <= NInput; i++)
  {
    toReturn += solution[i]*solution[i];
  }

  double toReturnS = sqrt(toReturn);
  return toReturnS;
}